

#Getting the local path to R libraries
if(Sys.info()[["sysname"]]=="Windows"){
	library(Biobase)
}else{
	locFile =file("../Misc/localization_files.txt","r")
	for(line in readLines(locFile)){
		if(substr(line,0,9) == "Rlib_path"){
			Rlib_path = gsub(" ","",sub("=","",sub("Rlib_path","",line)))
		}
	}
	library(Biobase, lib.loc=Rlib_path)
}



pngFile<-paste("../Images/genes_vs_genotypes.png",sep="")
expressionFile="../TXT/explvl.txt"
genotypeFile="../TXT/genotype.txt"
annotationFile="../TXT/annotation.txt"
tissue_nameFile="../TXT/tissuetype.txt"



#REFORMING FROM PHP
#loading data
require(Biobase) 
expression<-read.table(expressionFile,sep="\t",row.names=1,header=TRUE)
colnames(expression)<-sub("^X","",colnames(expression)) 

genotype<-read.table(genotypeFile,sep="\t",row.names=1,header=TRUE)
colnames(genotype)<-sub("^X","",colnames(genotype)) 

tissue_name<-read.table(tissue_nameFile,stringsAsFactors=FALSE)[,1]

annotation<-read.table(annotationFile,sep="\t",row.names=1,header=TRUE)


if(!all(rownames(genotype)==rownames(expression)))stop("rownames of genotypeFile and rownames of expressionFile must be equal")
if(!all(rownames(annotation)%in%colnames(expression)))stop("Some annotation rownames were not found in expression file")
if(!all(colnames(expression)%in%rownames(annotation)))stop("Some expression rownames were not found in annotation file")

#forming ExpressionSet
expressionset<- new(
		"ExpressionSet", 
		exprs = t(as.matrix(expression)), 
		experimentData = new("MIAME",title=tissue_name))
pData(expressionset)<-as.data.frame(genotype)






#doing statistics
garbage<-Sys.setlocale("LC_COLLATE", "C")
correlations<-data.frame(pvalue=vector(),estimate=vector())


for(snp in colnames(pData(expressionset))){
	for(probeset in featureNames(expressionset)){
		expression <- exprs(expressionset)[probeset,]
		genotype <- as.numeric(expressionset[[snp]])
		model <- lm(expression ~ genotype)
		correlations[paste(snp,probeset),"pvalue"] <- summary(model)[["coefficients"]]["genotype","Pr(>|t|)"]
		correlations[paste(snp,probeset),"estimate"] <- summary(model)[["coefficients"]]["genotype","Estimate"]
	}
}
correlations<-correlations[order(correlations[,"pvalue"]),]


#CREATING PNG
if(nrow(correlations)> 15){
	correlationsForPNG<-correlations[1:15,]
	numberOfRows <- 1 + ceiling(nrow(correlationsForPNG) / 3)
}else{
	correlationsForPNG<-correlations
	numberOfRows <- ceiling(nrow(correlationsForPNG) / 3)
}


png(filename = "../Images/genes_vs_genotype.png", width = 1440, height = 480*numberOfRows)
par(omi=c(0.2,0.2,0.2,0.2),mai=c(1,1,1,1))
layout(matrix(1:(3*numberOfRows),ncol=3,byrow=TRUE) )
for(i in 1:nrow(correlationsForPNG)){
	rowEntry<-rownames(correlationsForPNG)[i]
	snp<-sub(" .+$","",rowEntry)
	probeset<-sub("^.+ ","",rowEntry)
	
	plot(
			y=exprs(expressionset)[probeset,],
			x=as.factor(pData(expressionset)[,snp]),
			ylab=list(paste(probeset,": ",annotation[probeset,"genesymbol"],sep=""),cex=1.5),
			xlab=list(snp,cex=1.5),
			main=list(paste("Expression plot of",probeset,"varying with",snp),cex=1.5),
			sub=list(paste("pvalue:",signif(correlationsForPNG[rowEntry,"pvalue"],3),sep=" "),cex=1.5)
	)
	mtext(paste("       on Y-axis: ",annotation[probeset,"genesymbol"],sep=""),adj=0,padj=-1.4,cex=1.2)
	mtext(paste("       on X-axis: ",snp,sep=""),adj=0,padj=-0.2,cex=1.2)
}
if(nrow(correlations)> 15){
	frame()
	frame()
	mtext("more plots available in pdf file",cex=2)
}
garbage <- dev.off()



#CREATING PDF
pdf("../PDF/genes_vs_genotype.pdf", width = 3*5, height = 5*5)
layout(matrix(1:15,ncol=3,byrow=TRUE) )
par(omi=c(0.2,0.2,0.2,0.2),mai=c(1,1,1,1))
for(i in 1:nrow(correlations)){
	rowEntry<-rownames(correlations)[i]
	snp<-sub(" .+$","",rowEntry)
	probeset<-sub("^.+ ","",rowEntry)
	plot(		
			x=as.factor(pData(expressionset)[,snp]),
			y=exprs(expressionset)[probeset,],
			ylab=list(paste(probeset,": ",annotation[probeset,"genesymbol"],sep=""),cex=1.5),
			xlab=list(snp,cex=1.5),
			main=list(paste("Expression plot of",probeset,"varying with",snp),cex=1.5),
			sub=list(paste("pvalue:",signif(correlationsForPNG[rowEntry,"pvalue"],3),sep=" "),cex=1.5)
	)
	mtext(paste("       on Y-axis: ",annotation[probeset,"genesymbol"]," - ",substr(annotation[probeset,"genename"],1,40),sep=""),adj=0,padj=-1.4,cex=0.7)
	mtext(paste("       on X-axis: ",snp,sep=""),adj=0,padj=-0.2,cex=0.7)
}
garbage <- dev.off()










#CREATING TXT
pdata<-cbind(data.frame(genesymbol=rep("",ncol(pData(expressionset))),genename=rep("",ncol(pData(expressionset)))), t(pData(expressionset)),stringsAsFactors=FALSE)
expression<-cbind(annotation[featureNames(expressionset),],as.data.frame(exprs(expressionset)))
rawData<-rbind(
		pdata,
		expression
)
colnames(rawData)[1:2]<-c("","")
write.table(rawData, file="../TXT/raw_data.xls", sep="\t",col.names=NA)
